Asymmetry of the work probability distribution 
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We show, both analytically and numerically, that for a nonlinear system making a transition from 
one equilibrium state to another under the action of an external time dependent force, the work 
probability distribution is in general asymmetric. 

PACS numbers: 

We consider a system in contact with a heat bath, which is driven out of equilibrium by an external time dependent 
force. This force drives it from an equilibrium state A to another equilibrium state B. It was shown by Jarzynski 
P, 0, Q that the equilibrium free energy difference, AF between these states can be related to the probability 
distribution of the work done W in taking the system from A to B. In particular, 

-AF / -W 



e kt — \ e KT j (1) 

where K is Boltzmann constant and T is the temperature. A related equality is that due to Crooks [3, Q which 
relates the probabilities for the forward process (^4 — > B) to the backward process (B — > A). A couple of years ago a 
simple but effective experiment on a mechanical oscillator was done by Douarche, Ciliberto, Patrosyan and Rabbiosi 
(DCPR)[y, 0] who showed that for a Gaussian distribution of W, Eq (1) can be recast as 

AF = (W) - ^ K —1L± (2) 

This relation which is similar to a relation found by Landau and Lifshitz 8] in the context of thermodynamic 
perturbation theory, was demonstrated experimentally and also analytically for a particular kind of forcing by DCPR 
for linear oscillators. The particularly convenient form of Eq (2) made us investigate whether it holds for nonlinear 
systems. In what follows, we have looked at the general system (although in highly viscous limit). Very recently for 
a quadratic V(x), the system has been explored by Douarche et alQ. 

dV 

mx + kx = -— + M(t) + f(t) (3) 
ox 

where M(t) is an externally applied time dependent force and f(t) is a random force that allows the system to be 
in equilibrium in absence of M(t). The system is supposed to be in equilibrium (state A) at t = and then we switch 
on F(t) for a time r, after which F(t) takes a constant value F(t). The system reaches to state B and equilibriates. 
In going from state A to state B, the work done is 

W = - [ M(t)x{t)dt (4) 



We are interested in the moments of W. We will work in the highly damped limit where the inertial term in Eq.(3) 
can be dropped. For a quadratic V(x) (linear system), we will prove Eq.(2) for an arbitrary M(t) and then go on to 
show that Eq.(2) needs to be modified for arbitrary V{x). The most significant finding is that, even for a symmetric 
V(x), the odd moments of AW = W— (W) are non- vanishing and hence the distribution of AW is asymmetric for 
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all non-quadratic V(x). We will show this analytically, using perturbation theory; suggest a generalization of Eq.(2) 
and numerically establish that the probability distribution P(AW) is indeed asymmetric in general. 
We begin with the linear harmonic oscillator under the action of a deterministic force M(t) and random force f(t). 
In the highly viscous limit, the system evolves according to 

x + Tx = M(t) + f(t) (5) 

where the random force f(t) has the correlation function 

(f(t)f(t')) = 2KT8(t-t') (6) 

We calculate the moments of the work done, from above dynamics. The solution for x(t) can be written down as 

x(t) = J G(t-t')[M(t') + f(t')]dt' (7) 

where G(t — t') is the causal Green function, 

G{t-t') = Q{t-t')e- T{t - t ' ) (8) 

Clearly, the average of the work is 

(W)=- J rftiM(ti) J 1 G{h - t 2 )M(t 2 )dt 2 (9) 
while, the deviation from the average is 

AW = W-(W) = - J dhM{h) J 1 G(t 1 -t 2 )f{t2)dt 2 (10) 
which has the mean square value of 

/T PT 
dt x M(t x ) \ dt 2 M(t 2 ) (11) 

J' 2 dt" J'' dt'G(h - t')G{t 2 - t")5{t' - t") 

where Eq.(6) has been used. Noting the identity, derived from the time translational invariance 

/ 2 G(t 2 - t")G(h - t")dt" = G{h - t 2 )/2T (12) 
Jo 



we arrived at 



((AW)") 



1 J 7 dhMih) J 



dhMih) / dt 2 M(t 2 )G(h-t 2 ) (13) 



2KT V 

Using the above and Eq.(9) 

(W) - = - f dhMih) J* dt 2 G(h t 2 )[M(t 2 ) + ^p±] (14) 

Integrating by parts the first term in the integral above and using G(0) = due to causality, we find 

m _«mr>__* (15) 

The free energy change is precisely this amount, and that establishes Eq.(2) for an arbitrary forcing M(t). 
We now consider the inclusion of a non linear term in the motion, which becomes 

x + Tx + Xx 3 = M(t) + f{t) (16) 
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The question to ask is whether the equality in Eq.(2) still holds? To investigate this, we specialize to the case of 
M(t) = M t as studied by DCPR and carry out a perturbative calculation to 0(A). 
We write 



X = X + Xxi + \ 2 x 2 + ... 
and substituting in Eq.(16) and equating the coefficients of equal powers of A on either side we get 

xo + Txo = M{t) + f(t) 

x\ + r^i = — x\ 

and so on. We can now expand W(t) according to Eq.(4) and (17), 

W = W + XW 1 + ... 

where 



From Eq.(18) we get 



W = - J M(t)x (t)dt 
Wx = - J M(t) Xl (t)dt 

x = [ G(t-t')[M(t') + f(t')}dt' 
Jo 

Xl = - I G(t - t')x 3 dt' 
Jo 



(17) 
(18) 

(19) 
(20) 



(21) 



We note that in the way we set it up,G(ti — t 2 ) is exactly the same G that we had for the linear problem. We have 
already calculated (Wo) (Eq.(9)) and now we concentrate on (W\). We write 



/T ft ft 

Mdt / dt'G(t - t') / dhG(t' - ti)M(ti) 



+ 3 J Mdt J dt'G(t-t') J dt-tdtidtz 
M(h) G(t' - h)G{t' - t 2 )G(t' - t s ) (f(t 2 )f(t 3 )) 
The second term in the r.h.s vanishes once we use Eq.(12) and causality. We are left with 

rt' 



(22) 



(W) = / Mdt / dt'G(t - t') 



dhG{t' -ti)M(ti) 



We now use M(t) — M t and carry out the the integration to arrive at 



2t 3 15r 2 16r 
+ 



4r r 2 2r 3 r 4 



(23) 



(24) 



, keeping the leading order terms, i.e. terms which increase with r. 

Let us now turn to the calculation of the variance ((W — (W)) 2 ). The perturbative calculation generates the following 
upto 0(A) 



((AW) 2 ) = ((Al^o) 2 ) + 2A (AWoAWi) 



(25) 



Where A Wo = Wq — (Wo) and AW\ = W\ — (W\). We have already calculated the first term in the r.h.s. Now we will 
concentrate on the second term. In the O(X) correction in the variance, we have found that the disconnected parts 
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(i.e. where the averaging is over AIT and AWi separately.) do not contribute. Specializing to the case M(t) = Mot, 
calculation leads to the 0(A) term of the variance 



Thus, corrected upto 0(A) we get 



0(A) part of 



((AW) 2 ) 
2KT 



Ml 

r 3 



2r3 

r 2 



27t 2 16t 

^ 3_ + T r 



, keeping the leading order terms. Thus, 



2KT 



2KT 



Hence following Eq.(2) and above, we see that (i) AF is no longer (W) - 



4r 4 

((AW) 2 
2KT 



r e 



(26) 



(27) 



(28) 



■, and (ii) for M = Mot, the difference 



arises at 0(t ) and 0(r ). (i) implies a non-Gaussian distribution for TT which perhaps not so surprising but it is 
(ii) which has the interesting consequence, (ii) leads to the conclusion that correction to the result has to have an 
asymmetric part. This is because, if the correction to the Gaussian distribution is symmetric, then the first correction 
that would arise is the flatness factor (((AIT) 4 ) — 3 ((AIT) 2 ) ), but a rather long and tedious calculation reveals that 
leading contribution to the term is at 0(t). The contribution at 0(r 2 ) can come only from ((ATT) 3 ). Does dynamics 
contribute to ((AIT) 3 ) ? We note that 



((ATT) 3 ) = 3A ((AIT ) 2 (AITi)) 



(29) 



and a cursory inspection of the solution of the equation of motion, we find that ((ATT) 3 ) is nonzero and the leading 
term is 0(t 2 ). 

This observation inspires us to start with a work probability distribution P(W) given by 



P(W) oc e [ 



(30) 



where [i 2 and a are parameters. Now we work out the expectation value of e kt with the above distribution 
keeping fi\, fi 2 small, and thus we have found 



exp 
1 + 



TT 



= e -[» + 



KT 

W (o I H 
KT\ {KT) 2 



2{KT) 2 ' x 



(31) 



M2 



(KT) 2 



(KTy 



, keeping only the linear order terms of [i\ and [i 2 . We now calculate different expectation values, like ((TT — TTo)), 
((TT — TTo) 2 ), ({W — TTo) 3 ) and ((TT — TTo) 4 ) and replacing m, ^ and a in above by these expectation values, we 
get the following basic result [using Eq.(2)] 



AF 



(TT) 



(TT — (TT)) 2 (TT — (TT)) 3 



2KT 



6(KT) 2 



(32) 



1 

+ 24 



3((TT-(TT)) 2 ) 2 -((TT-(TT)) 4 ) 



(KT) 3 

To test the relation numerically, we decided to work first with the quadratic nonlinearity in equation of motion, 

x + Tx + Xx 2 = M(t) + f(t) (33) 

We will restrict ourselves to such values of A with r = 1 (we take T = 1 every where), so that trajectory does not run 
away. In this case, the dynamics does not generate the last term in Eq.(32) to the lowest order in A. It is the cubic 
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deviation which is the most significant and that would imply an asymmetric probability distribution for the work W . 
This is not unexpected since the potential for Eq.(33) is cubic and hence asymmetric. We have taken, 

M{t) = t < (34) 
= M a t/r < t < t/M 
= 1 t > t/Mq 

We allow the system to equilibriate at t = and at t = by running it with M=0 for an interval and with M = 1 
for a similarly long interval. In between we generate the values of x at different points by the following, 

x(t + At) =x(t)- (x(t) + Xx 2 (t))At + M(t)At + V2KTAtr](t) (35) 

where r](t) is a random number between and 1. We calculate all the quantities in the unit of 2KT. We calculate 
a trajectory [a;(£)] T starting from an initial value and evaluate the work according to Eq.(4). The ensemble is one of 
initial conditions x(0) and we calculate (W), ((AM 7 ) 2 ), ((AW) 3 }. We have also found the work distribution function 
P(W). For A = (i.e. the system with quadratic potential) and for A = 20 the distributions are shown in fig. (1,3) 
respectively. For nonzero A the distribution is asymmetric, which is expected for an asymmetric potential. The free 
energy difference between times t = and t = r can be found as 

AF = -x 2 + -x 3 (36) 
2 3 K J 

It does not match with the r.h.s of Eq.(2) for nonzero A instead we calculate ((AW 7 ) 3 ) and find that Eq.(36) gives a 
correct description. 

We repeated the numerics with quartic oscillator (i.e.V(x) = \x 2 + jx 4 )a.nd found that ((AW) 3 ) is certainly nonzero, 
indicating an asymmetric distribution. For small values of A, the asymmetry is striking. For large values of A due to 
dominating ((AW) 4 ) the distribution becomes sharply peaked, and the asymmetry is difficult to make out, although 
its existence is guaranteed by the nonzero value of ((AW) 3 ). The distribution for A = 0.1 and A = 20 is shown in 
fig. (2,4) respectively. In fig(5) comparison between the work distribution functions for the cases when A = 0.1 and 
A = are shown by plotting together. The asymmetry is clear from it. Recently Mai and Dhar [ujhave found an 
asymmetric distribution for V(x) — ax 2 + bx 3 + cx 4 . Our contention is that the asymmetry exists even if b = 0. Here 
we have to calculate AF by 

AF = -x 2 + -x 4 (37) 
2 4 

We have calculated AF using Eq.(32) also, after calculating required moments from the work probability distribution 
obtained, and have found again, that it gives the correct description. 

We end by pointing out a possible application. We consider a ferromagnet or an Ising magnet near but above its 
critical point T c . We can imagine being close to T c , but sufficiently far away so that the mean field Landau model 
is valid. If we now switch on an time dependent magnetic field, then the dynamics of the mean magnetisation will 
be given by an equation of the form shown by Eq.(16). If we are in the region T < T c , then the dynamics will be 
governed by Eq.(16) with an added quadratic nonlinearity- the kind considering Ref.10. It will be interesting to check 
the veracity of Eq.(32) in this case. 
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FIG. 1: Here force is applied for 8.4 sec. Mo and r are 100 and 840 respectively, in S.I. units. Here A = 0. Note the symmetry 
of the distribution 
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FIG. 2: Here force is applied for 8.4 sec. Mo and r are 100 and 840 respectively, in S.I. units. Here A = 0.1. The tail in left 
side signifies the asymmetrical nature of the distribution here. Here V[x) = (l/2):r 2 + (0.1/4)x 4 . 
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FIG. 3: Here force is applied for 8.4 sec. Mo and r are 100 and 840 respectively, in S.I. units. Here A = 20. The tail in left 
side signifies the asymmetrical nature of the distribution here. Here V(x) = (l/2)x 2 + (20/3)a; 3 . 




FIG. 4: Here force is applied for 8.4 sec. Mo and r are 100 and 840 respectively, in S.I. units. Here A = 20. The symmetrical 
nature is due to dominance of ((W — (W0) 4 ) °f tne distribution here. Here V(x) = (l/2)x 2 + (20/4)a; 4 . 
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FIG. 5: For same values of Mo, r and the time of application of M(t) as before, but for V(x) = (l/2)x 2 and V(x) = 
(l/2)x 2 + (0.1/4)x 4 the above liny and dotted plots are done respectively. The striking asymmetry in nonlinear case is clear 
here. 



